function dataTrend = TrendExt(res,unm)
%TrendExt - Hector/CATS-based trend evaluation and data extraction
%
%       dataTrend = TrendExt(RES,UNM)
%
% Let RES be the result of a system session aimed at estimating the 
% trend of a component of a GNSS time series by means of Hector 
% (for example, if in Windows the Hector directory is C:\GNSS\hector
% and the control file is estimatetrend.ctl, placed in the same directory 
% of MATLAB, the command is 
%   [STATUS,RES] = system(...
%   'wsl /mnt/c/GNSS/hector/estimatetrend estimatetrend.ctl'),
% where the fact that Hector operates in LINUX is taken into account).
% This function extracts the trend data providing the struct variable 
% dataTrend from which the data can be easily transferred to a tsData 
% object.
% If UNM is undefined, empty, or is a string different from 'mm' or 'm',
% the unit is millimeter (string 'mm'). If UNM is defined and is either
% 'mm' or 'm', such a unit is used.

% G. Teza, 2021

if (nargin < 2) || isempty(unm) || ~ismember(unm,{'mm','m'})
    unm = 'mm';
end

dataTrend = struct('unit',unm); 

dataTrend.trend = [];
dataTrend.strend = [];
dataTrend.bias = [];
dataTrend.sbias = [];
dataTrend.datebias = [];

dataTrend.wnfraction = [];
dataTrend.wnsigma = [];

dataTrend.plfraction = [];
dataTrend.plsigma = []; 
dataTrend.plk = []; 
dataTrend.plsk = []; 

dataTrend.cosy = [];
dataTrend.scosy = [];
dataTrend.siny = [];
dataTrend.ssiny = [];
dataTrend.ampy = [];
dataTrend.sampy = [];
dataTrend.phay = [];
dataTrend.coshy = [];
dataTrend.scoshy = [];
dataTrend.sinhy = [];
dataTrend.ssinhy = [];
dataTrend.amphy = [];
dataTrend.samphy = [];
dataTrend.phahy = [];

It1 = strfind(res,'bias');
if ~isempty(It1)
    It2 = strfind(res,'trend');
    if numel(It2) > 1
        [~,nt2] = min(abs(It2-It1));    % search of the correct index
        It2 = It2(nt2);
    end
    It3 = strfind(res,[unm '/year']);
    % bias:
    resb = res(It1:It2-2);
    Jb1 = strfind(resb,':');
    Jb2 = strfind(resb,'+/-');
    Jb3 = strfind(resb,unm);
    Jb4 = strfind(resb,'at');
    Jb5 = strfind(resb,',');
    bias = str2double(resb(Jb1+1:Jb2-1));
    biass = str2double(resb(Jb2+3:Jb3-1));
    dates = resb(Jb4+3:Jb5-1);
    ys = str2double(dates(1:4));
    if strcmpi(dates(8),'/')
        ms = str2double(dates(6:7));
        ds = str2double(dates(9:end));
    else
        ms = str2double(dates(6));
        ds = str2double(dates(8:end));
    end
    dn = datenum(ys,ms,ds);
    dateb = serial2frac(dn);
    % trend:
    rest = res(It2:It3-1);
    Jt1 = strfind(rest,':');
    Jt2 = strfind(rest,'+/-');
    trend = str2double(rest(Jt1+1:Jt2-1));
    trends = str2double(rest(Jt2+3:end));
    dataTrend.bias = bias;
    dataTrend.sbias = biass;
    dataTrend.datebias = dateb;
    dataTrend.trend = trend;
    dataTrend.strend = trends;
else
    warning('No valid trend');
    return
end

Ipl1 = strfind(res,'Powerlaw:');
if ~isempty(Ipl1)
    Ipl2 = Ipl1+150;
    respl = res(Ipl1:Ipl2);
    Jpl1 = strfind(respl,'fraction');
    Jpl2 = strfind(respl,'sigma');
    Jpl3 = strfind(respl,'d  ');
    Jpl4 = strfind(respl,'kappa');
    resplf = respl(Jpl1:Jpl2-1);
    Kplf1 = strfind(resplf,'=');
    fraction = str2double(resplf(Kplf1+1:end));
    respls = respl(Jpl2:Jpl3-1);
    Kpls1 = strfind(respls,'=');
    Kpls2 = strfind(respls,unm);
    sigma = str2double(respls(Kpls1+1:Kpls2-1));
    respld = respl(Jpl3:Jpl4-1);
    Kpld1 = strfind(respld,'+/-');
    ds = str2double(respld(Kpld1+3:end));
    kappas = 2*ds;
    resplk = respl(Jpl4:end);
    Kplk1 = strfind(resplk,'=');
    Kplk2 = strfind(resplk,'+/-');
    kappa = str2double(resplk(Kplk1+1:Kplk2-1));
    dataTrend.plfraction = fraction;
    dataTrend.plsigma = sigma;
    dataTrend.plk = kappa;
    dataTrend.plsk = kappas;
end

Iwn1 = strfind(res,'White:');
if ~isempty(Iwn1)
    Iwn2 = Iwn1+60;
    reswn = res(Iwn1:Iwn2);
    Jwn1 = strfind(reswn,'fraction');
    Jwn2 = strfind(reswn,'sigma');
    Jwn3 = strfind(reswn,unm);
    reswnf = reswn(Jwn1:Jwn2);
    Kwnf = strfind(reswnf,'=');
    wnf = str2double(reswnf(Kwnf+1:end-1));
    reswns = reswn(Jwn2:Jwn3);
    Kwns = strfind(reswns,'=');
    wns = str2double(reswns(Kwns+1:end-1));
    dataTrend.wnfraction = wnf;
    dataTrend.wnsigma = wns; 
end

dataTrend = OscillExt(res,dataTrend,'yearly',unm);  % ancillary function, see below
dataTrend = OscillExt(res,dataTrend,'hyearly',unm);


%% ANCILLARY FUNCTIONS:

function dataTrend = OscillExt(res,dataTrend,str_y,unm)
% Yearly/hyearly oscillation data extraction
Iy1 = strfind(res,['cos ' str_y]);
if ~isempty(Iy1)
    Iy2 = strfind(res,['sin ' str_y]);
    Iy3 = strfind(res,['Amp ' str_y]);
    Iy4 = strfind(res,['Pha ' str_y]);
    Iy5 = strfind(res,'degrees');
    if numel(Iy5) > 1
        [~,ny5] = min(abs(Iy5-Iy4));    % search of the correct index
        Iy5 = Iy5(ny5);
    end
    % cos:
    rescy = res(Iy1:Iy2-2);
    Jcy1 = strfind(rescy,':');
    Jcy2 = strfind(rescy,'+/-');
    Jcy3 = strfind(rescy,unm);
    cosyv = str2double(rescy(Jcy1+1:Jcy2-1));
    cosys = str2double(rescy(Jcy2+3:Jcy3-1));
    % sin:
    ressy = res(Iy2:Iy3-2);
    Jsy1 = strfind(ressy,':');
    Jsy2 = strfind(ressy,'+/-');
    Jsy3 = strfind(ressy,unm);
    sinyv = str2double(ressy(Jsy1+1:Jsy2-1));
    sinys = str2double(ressy(Jsy2+3:Jsy3-1));
    % amplitude:
    resay = res(Iy3:Iy4-2);
    Jay1 = strfind(resay,':');
    Jay2 = strfind(resay,'+/-');
    Jay3 = strfind(resay,unm);
    ampyv = str2double(resay(Jay1+1:Jay2-1));
    ampys = str2double(resay(Jay2+3:Jay3-1));
    % phase
    respy = res(Iy4:Iy5+6);
    Jpy1 = strfind(respy,':');
    Jpy2 = strfind(respy,'degrees');
    phayv = str2double(respy(Jpy1+1:Jpy2-1));
    if strcmpi(str_y,'yearly')
        dataTrend.cosy  = cosyv;
        dataTrend.scosy = cosys;
        dataTrend.siny  = sinyv;
        dataTrend.ssiny = sinys;
        dataTrend.ampy  = ampyv;
        dataTrend.sampy = ampys;
        dataTrend.phay  = phayv;
    elseif strcmpi(str_y,'hyearly')
        dataTrend.coshy  = cosyv;
        dataTrend.scoshy = cosys;
        dataTrend.sinhy  = sinyv;
        dataTrend.ssinhy = sinys;
        dataTrend.amphy  = ampyv;
        dataTrend.samphy = ampys;
        dataTrend.phahy  = phayv;
    end
end